#ifndef SCALEFACTOR_H
#define SCALEFACTOR_H

//#include <iostream>
#include <stdio.h>
#include <math.h>
/*
class TN {
		public:
			int N[7];
			int &operator =(int i)   {  N[0]=i; return i;  }
			int &operator [](int i)  {  return N[i+1];     }
			operator int(){	return N[0]; }
} N;
*/

class Value{
	public:
		double value;
		double &operator =(double value)   {  this->value = value; return value;  }
		operator double(){	return value; }
		Value (): value(0) {};
		Value (double v): value(v) {};

};

class Sec:public Value{
	public:
		
		using Value::operator =;
		using Value::operator double;
		Sec(): Value(){}
		Sec(double v): Value(v){}
		
		friend Sec operator -(Sec sec1, Sec sec2){
			return sec1.value-sec2.value;
		};

		double Gyr(){ return  value*3.16880e-17;  }
		double Myr(){ return  value*3.15576e+13;  }
};





class ScaleFactor: public Value {
	public:
		
		using Value::operator =;
		using Value::operator double;
		
		ScaleFactor(): Value(){}
		ScaleFactor(double v): Value(v){}
		
		
		double z(){
			
			
		}
		
		void z(double z){
			value = 1/(z+1);
		}
		
		Sec s(double Omega0, double OmegaLambda, double h){
			
			if(OmegaLambda+ Omega0!=1) { printf("only implemented for flat cosmology so far.\n"); return -1; }
			
			double H0 = 3.2407789e-18; // h/sec
			double a = value;
			
			double factor1 = 2.0 / (3.0 * sqrt(OmegaLambda));

			double term1 = sqrt(OmegaLambda / Omega0) * pow(a, 1.5);
			double term2 = sqrt(1 + OmegaLambda / Omega0 * pow(a, 3));
			double factor2 = log(term1 + term2);

			return (factor1*factor2)/(H0*h);
			
		}
		

		Sec H(double Omega0, double OmegaLambda){
			
			if(OmegaLambda + Omega0 != 1) { printf("only implemented for flat cosmology so far.\n"); return -1; }
			
			double H0 = 3.2407789e-18; // h/sec
			double a = value;
			return H0*sqrt(Omega0*a + OmegaLambda*pow(a,4))/pow(a,2); // only flat case
			
		}

		
/*		
		
		
		double H0;
		double HubbleParam;
		double Omega0;
		double OmegaLambda;
		double SecPerMegayear;
		double pc;
		
		TimeConverter():
		H0(3.2407789e-18),  // in h/sec 
		HubbleParam(0.73),
		Omega0(0.25),
		OmegaLambda(0.75),
		pc(3.08568025e16), // in meters (SI)
		SecPerMegayear(3.155e13)
		{}
		
		//  in seconds

		
		double da2dt(double a1, double a2){
			return a2t(a2)-a2t(a1);
		}
		
		double H(double a){
			return H0*sqrt(Omega0*a + OmegaLambda*pow(a,4))/pow(a,2); // only flat case
		}
		*/
};



/*

double a2t(double a, double Omega0, double OmegaLambda, double HubbleParam){
			
	if(OmegaLambda + Omega0 != 1) { printf("only implemented for flat cosmology so far.\n"); return 0; }
		
	double factor1 = 2.0 / (3.0 * sqrt(OmegaLambda));
			
	double term1 = sqrt(OmegaLambda / Omega0) * pow(a, 1.5);
	double term2 = sqrt(1 + OmegaLambda / Omega0 * pow(a, 3));
	double factor2 = log(term1 + term2);

	return (factor1*factor2)/(H0*HubbleParam);
		
}


double a2t(double a, double Omega0, double OmegaLambda, double HubbleParam){
			
	if(OmegaLambda + Omega0 != 1) { printf("only implemented for flat cosmology so far.\n"); return 0; }
		
	double factor1 = 2.0 / (3.0 * sqrt(OmegaLambda));
			
	double term1 = sqrt(OmegaLambda / Omega0) * pow(a, 1.5);
	double term2 = sqrt(1 + OmegaLambda / Omega0 * pow(a, 3));
	double factor2 = log(term1 + term2);

	return (factor1*factor2)/(H0*HubbleParam);
		
}


*/








#endif